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Abstract 

The linear response description for impurity diffusion in a granular fluid undergoing homoge- 
neous cooling is developed in the preceeding paper. The formally exact Einstein and Green-Kubo 
expressions for the self-diffusion coefficient are evaluated there from an approximation to the ve- 
locity autocorrelation function. These results are compared here to those from molecular dynamics 
simulations over a wide range of density and inelasticity, for the particular case of self-diffusion. It is 
found that the approximate theory is in good agreement with simulation data up to moderate den- 
sities and degrees of inelasticity. At higher density, the effects of inelasticity are stronger, leading 
to a significant enhancement of the diffusion coefficient over its value for elastic collisions. Possible 
explanations associated with an unstable long wavelength shear mode are explored, including the 
effects of strong fluctuations and mode coupling. 



I. INTRODUCTION 



Attempts to describe granular media in terms of a more fundamental underlying statistical 
mechanics have met with considerable success. As a prototype for this approach, in the 
preceding paper standard linear response methods from normal fluids have been applied, 
mutatis mutandis, to the case of an impurity particle diffusing in an isolated one component 
fluid of smooth inelastic hard spheres {d = 3) or disks {d = 2) of diameter a. The collisions 
are characterized by a coefficient of normal restitution a. The isolated system (or with 
periodic boundary conditions), is not in the typical Gibbs state as for elastic collisions, but 
rather in a time dependent homogeneous cooling state (HCS). It has been shown in and 
will be elaborated further here, that this time dependent HCS can be exactly transformed 
to a stationary state description. In the present work, this stationary state description is 
evaluated by molecular dynamics (MD) simulation to measure the mean square displacement 
of the impurity, its velocity autocorrelation function, and the resulting diffusion coefficient 
defined in terms of a formal Einstein or Green-Kubo relation, respectively. For practical 
purposes, attention is restricted to the case of self-diffusion, where the mechanical properties 
of the impurity are the same as those of the fluid particles. 

The velocity autocorrelation function was approximated in using a cumulant expansion 
and also by means of kinetic theory methods. In the usual first Sonine polynomial expansion 
for the latter, this leads to a simple exponential decay in an appropriate dimensionless time 
scale defined below. For the numerical results considered here, one further approximation is 
made in evaluating the decay rate, namely two-particle velocity correlations are neglected. 
In this case, the results are the same as those obtained from the Enskog kinetic theory. 
Earlier studies of self-diffusion show excellent agreement of this theory with MD results 
at very low number density n, over a wide range of values for a. The present work extends 
that study to higher densities. More specifically, three-dimensional systems with densities 
in the interval 0.1 < n* = na^ < 0.75, for 0.5 < a < 1, will be considered. 

The HCS is known to be unstable under long wavelength perturbations or fluctuations . 
To avoid this problem, the system size in most cases considered is chosen to be smaller than 
the critical wavelength. The latter is a function of the density and inelasticity, decreasing 
with increasing values of each. Consequently, at the highest densities and smallest a a small 
system of 108 particles is required. At moderate densities and inelasticities, a system size 
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of 512 particles is used. The agreement of theory and simulation is found to be quite good 
for all a at n* < 0.25 and for all densities at a > 0.9. At higher densities and smaller a, 
significant discrepancies occur, with the diffusion coefficient obtained from MD being almost 
an order of magnitude greater than the theoretical estimate at n* = 0.75 and a = 0.5. The 
failure of the Enskog theory at high densities is well-known for the case of elastic collisions, 
due to cage effects and correlated binary collisions. However, in that case the diffusion 
coefficient from MD is smaller than the theoretical prediction. Therefore, it is clear that a 
quite different density mechanism is effective for inelastic collisions. 

It is reasonable to expect that the underlying instability is responsible in part for the 
above discrepancies. Although the instability is, in principle, avoided by the use of small 
systems, the latter is prone to large fluctuations. These are not quantitatively significant at 
a = 1 but grow rapidly with decreasing a, as can be seen in the noise level for the kinetic 
temperature. This will be illustrated later on in Fig. |I[ A related issue is the role of the 
instability in affecting the usual mode coupling corrections to the Enskog theory. Such mode 
coupling terms arise from correlations of spontaneous fluctuations in the long wavelength 
hydrodynamic modes, which are significantly modified by the potential instability. Enhanced 
fluctuations will be addressed in the qualitative analysis of the data below. 

The plan of the paper is as follows. In Sec. II, the relevant results in the preceding paper 
are shortly summarized, and a new scale transformation, useful for practical purposes, is 
introduced. The equations of motion of the new phase-space variables lead in a direct way to 
a steady-state simulation method for the HCS. Starting from an arbitrary initial condition, 
the system rapidly approaches a steady state, whose properties are simply related to those 
of the HCS. It is important to stress that this relationship is not approximate, but an 
exact consequence of a change of variables. In particular, the method to measure the mean 
square displacement of a tagged particle as well as its velocity autocorrelation function is 
discussed in detail. Moreover, both quantities are shown to lead to equivalent results for the 
self- diffusion coefficients. Although in practice some discrepancies appear in the numerical 
results obtained by the two procedures, their origin is well understood. 

Comparison of theory and simulation is also started in Sec. II, and completed in Sec. |T|. 
As mentioned above, the agreement is fairly good at low densities, but relevant effects, which 
are not taken into account by the theory developed in [|^], show up even at moderate density. 
The peculiar nature of these effects for inelastic collisions is discussed. Numerical evidence 
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is provided indicating that the underlying hydrodynamic instabihty associated to the shear 
mode plays an essential role in the "anomalous" behavior of the self-diffusion coefficient. 
Nevertheless, such changes in the diffusion constant do not compromise the existence of the 
diffusion process, which is confirmed for such conditions by simulation results for the mean 
square displacement. Finally, Sec. |V] contains a short summary of the results and also some 
indications of the possible extensions of the reported work. 

II. STEADY-STATE SIMULATION METHOD 

To investigate the nature of diffusion in the HCS and to test the theoretical results 
presented in the preceding paper a more fundamental description via MD simulation 
will be considered. The HCS distribution function for a fluid of + 1 particles has the 
scaling form 

Here and in the following the notation is the same as in ref. ||T||. In particular, v{t) is the 
thermal velocity defined in the usual way (with Boltzmann's constant set equal to one) 
and i is proportional to the mean free path of the gas. The time dependence of the above 
distribution is due to coUisional cooling and is determined from 

dMt) = -Icmt), (2) 

where ({t) is the cooling rate. A direct simulation of the cooling fluid, as described by the 
Liouville dynamics in the actual phase variables, is difficult, since the rapid cooling of the 
fluid leads to numerical inaccuracies very soon. One method of dealing with this would be 
to periodically redefine the time-scale of the simulation, so that the typical particle velocity 
remains of the order of unity. However, if the rescaling is only used to control the numerical 
stability of the simulation, the state being simulated would nevertheless be time- dependent. 
For this reason, some simulation studies employ different type of thermostats, such as an 
externally imposed Brownian force or thermal boundaries, in order to generate a steady 
state. While these methods provide more or less realistic models of various experimental 
procedures, they obviously probe a state which is in some way related to, but not identical 
with, the homogeneous cooling state. 
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Soto et al. 0], noting the fact that there is no intrinsic time scale in the dissipative 
hard sphere model, have proposed to rescale all particle velocities after every collision, thus 
establishing a steady state similar to that described in [Q. However, this procedure has the 
effect of replacing the binary collision dynamics of the dissipative hard sphere model by an 
A^-body dynamics, since the result of a collision is to alter all atomic velocities and not just 
those of the colliding atoms. It is reasonable to imagine that the binary dynamics is recovered 
in the infinite system limit, but the connection between the HCS and this dynamics is not 
clear for the case of the small systems considered in most of the MD simulations. Instead, 
we follow the procedure used in ref . |^ , according to which an exact mapping of the HCS 
onto a steady state is exploited as the basis of the simulation method. Following the ideas 
developed in [|l|, all velocities are scaled relative to v{t) and the new dimensionless time is 
given by 

ds{t) = v{t)dt/i. (3) 

This time scale is a measure of the average collision number. The corresponding Liouville 
equation in these variables supports a stationary HCS solution given by p^^^ ({ct*j,v*}), 
where q* = q^/f and v* = \i/v{t). Moreover, the average value of any phase function 
A({q, v}) is given by Eq.(41) of reference [Q. 

{A-t) = I drp*{r)A{{eci* is),vm (.)}). (4) 

The dynamics in the new phase space is obtained from 

dsci: (s) = V* (.) , S.v* (.) = lc*v* (s) + L*v* (.) , (5) 

where L* is the dimensionless Liouville operator, and (* = i({t)/v{t) is the dimensionless 
cooling rate. This is the same as the usual dynamics for hard spheres, except that the term 
proportional to (*/2 represents an acceleration between collisions which balances the energy 
lost during collisions, thus enabling a steady state. 

Although it is possible to relate average values of the relevant dynamical functions of the 
original variables with the average values of the same functions of the scaled variables for a 
general situation, attention will be restricted in the following to the HCS. For homogeneous 
systems, (* becomes time independent. In principle, the scaled dynamics defined by Eq. 
(I) can be simulated for the properties of interest, without the complications of continuous 
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cooling in real time. In this formulation, the simulation is qualitatively similar to that for 
elastic collisions, in the sense that the system rapidly approaches a state for which the time- 
average of the instantaneous scaled temperature is constant, and subsequent determination 
of average properties is simpler and numerically more accurate. A technical complication is 
the need to know the exact value of (* a priori, which in general is not possible, since it is 
determined by the original dynamics. Consequently, it is useful to make a second change of 
scale, 

q- = q*, vr = K/av*, s* = {C/w*)s, A** (s*) = A (s) , (6) 

where w* is an arbitrary time independent dimensionless frequency. Equation (|^) then 
becomes 

with L** the same L*, but with the replacements {qi,v*} {^i*,^**}- This is the form 
most convenient for MD simulation, with any reasonable choice for w* . More explicitly, 
the particle dynamics implied by Eq.(|^) and implemented in our event-driven simulations 
consists of an accelerating streaming between collisions 

9,*qr = vr, (8) 

while the effect of the collision of two particles is to alter their relative velocity according to 

g^^ - = g^T - (1 + «) (q^ ■ g^) q^, (9) 

where g** = v** — v**, q** is the unit vector pointing from atom i to atom j and the center 
of mass velocity remains the same. 

The relationship of C* to w* is determined by the steady state temperature obtained in the 
simulation. To see this note that while the instantaneous kinetic energy of the entire system 
is clearly not constant, a corresponding "temperature" defined as T** = {vl*"^] s*)** does 
approach a constant, as is shown in Appendix ^, 



Thus (* = w* / a/2T**(oo) which relates quantities measured during the steady state simu- 
lation to those calculated in In summary, simulation of Eqs.(^ is expected to yield a 
steady state after a short transient period. Subsequently, ensemble averages of properties 
can be determined as time averages by making the usual assumption of ergodicity. These 
properties are directly related to those of the HCS by a simple scale transformation as 
described above. 

This steady-state simulation method removes any limitation on the time for which tra- 
jectories may be followed, but there are other limitations due to a long wavelength hydrody- 
namic instability]^] for systems with dimension larger than a critical size Lc = 27ri^y2ri* /(* 
where r]* = r]{t)/n'mi'{t)i is the shear viscosity. These instabilities occur when the decay rate 
of a shear mode fluctuation is less than the cooling rate of the thermal velocity so that their 
size grows relative to the temperature. In the scaled dynamics, the instabilities show up as 
ordinary unstable modes that grow exponentially with s*, and are therefore easily identified 
in the simulations 0. Since these growing modes represent a spontaneous breaking of the 
(assumed) spatial homogeneity of the system, the impurity particle no longer undergoes sim- 
ple diffusion and, therefore, its motion is beyond the scope of this paper. The critical size Lc 
is a function of the density and coefficient of restitution, being smaller for high density and 
small restitution coefficient. In all the cases reported here, the system size is smaller than 
that required for the instability to show up. For a system of 108 atoms and n* < 0.5, this 
provides no limitation on the accessible densities and coefficients of restitution, otherwise 
only sufficiently low (high) values of the density (coefficient of restitution) can be studied. 
A different limitation on the simulation is the "inelastic collapse" [^], which occurs for small 
values of a thereby setting a lower limit on the values accessible to simulation. For the 
system sizes considered in this paper, the above mentioned limits have been characterized 
in more detail elsewhere 0. 

The three-dimensional simulations reported in the following begin in all cases with an 
equilibrated fluid subject to elastic colhsions (i.e. a = l,w* = 0) with T** = 1/2. The 
coefficient of restitution is then set to the desired value and the scaling parameter w* is set 
to w* = CI; («) a/2T**(0) using the Enskog estimate for the cooling rate ("^ (a) as given by 
Eq. (74) in ref. ||l|. This choice ensures that the temperature remains of order 1 for all values 
of a, for optimal numerical stability, as well as provides a continuous path to the equilibrium 
system for which a = 1 and w* = 0. The simulation is then continued for 10^ collisions to 
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allow the system to reach the steady state at T** = w*'^/2(*'^ where (* is the true cooling 
rate generated by the simulation. Since the largest system considered consists of 500 atoms, 
this corresponds to at least 2 x 10'^ collisions per atom, and is more than sufficient to reach 
the steady state. A final simulation of 10'' additional collisions is then performed, during 
which all statistical averages of interest are accumulated. Figure shows the behavior of 
the instantaneous kinetic energy (or temperature) during several typical simulations. In all 
cases, the instantaneous temperature fiuctuates around a stationary average value, with the 
size of the fiuctuations increasing as the coefficient of restitution decreases. The fact that the 
average value itself varies with a is due to (a) the inadequacy of the simple estimate of (* we 
have used for decreasing values of a and (b) to the presence of long-lived shear fiuctuations, 
or vortices, which decay on a time scale of thousands of collisions near the instability. As 
will be discussed below, the latter give rise to large fiuctuations in all components of the 
kinetic contribution to the pressure tensor as well as raising the apparent temperature. 

The system is expected to show a typical diffusive behavior under the scaled dynamics 
after a short transient time, i.e. for s* » 1. Then the diffusion coefficient can be obtained 
either from the mean squared displacement (msd) using the Einstein relation, Eq. (70) of 
ref. or from the Green-Kubo relation, Eq. (69) of ref. |Q. The former is evaluated, for 
d = 3, from 

^ ^ ^ y \cC*{s*) - qr(0)|' . (11) 

Figure ^ shows some typical simulation data exhibiting the expected linearity of the mean- 
square displacement as a function of s*. The Green-Kubo relation expresses the diffusion 
coefficient in terms of the velocity autocorrelation function (vacf) as 

D*\s*) = - [ ds*' (v**(s*') • V**)**. (12) 
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The vacf is evaluated from 

N+l imax 



(JA), (13) 



where a, (3 = x, y, z, A is the sampling period, and jmax is a function of the amount of data 
available and the time at which the steady state is established. Notice that the stationarity 
of the scaled dynamics is explicitly used to evaluate the vacf for a given time separation s*, 
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by averaging over many different samplings during the simulation. In the simulations, A is 
taken to be 1/4 of the Boltzmann mean free time at T** = 1/2, i.e. A** = ^ {Ana^y/ny^ . 
As discussed above, T** = 1/2 is not the exact value of the converged temperature, but in 
formulating A** only a reasonable order of magnitude is required. To evaluate the right side 
of Eq. (|13D, we follow the procedure of Futrelle and McGinty [0] (see also the discussion in ref. 
0), according to which the components of the velocity for each atom sampled at intervals of 
A** are stored. When the number of values stored reaches a threshold, N^, the convolution 
theorem is used to evaluate Eq. (|T3D by means of the fast Fourier transform, and so to obtain 
an estimate of the velocity autocorrelation function. This procedure is then repeated until 
the end of the simulation and the various estimates, typically between tens and hundreds of 
samples, are averaged to obtain the final estimate of the correlation function. Figure ^ shows 
some typical correlation functions obtained in this way, using A^a = 1024. The functions 
have been normalized in all cases by their initial values, i.e. the function C**, as defined 
in the preceding paper, has been plotted. To get the diffusion constant, the integral of the 
velocity correlation function is computed using Simpson's rule and the variance of the pool 
of estimates is used as the basis for the calculation of the standard error. 

Both methods, Einstein relation and Green-Kubo relation, should yield the same values 
of the coefficient of self-diffusion. Figures 2 and 3 show that the mean square displacement 
becomes linear in s* and the velocity autocorrelation function decays to zero for s* ^ 10. 
For larger s* both forms for the diffusion coefficient approach a constant, and all further 
discussion is restricted to this diffusion constant. The consistency between the Green-Kubo 
and Einstein relations is confirmed in Fig. ^ where the estimates obtained by both proce- 
dures are compared for a system of 108 atoms and several values of the density and of the 
coefficient of restitution. While in qualitative agreement, the difference in values obtained 
by the two methods is in some cases as great as 10%. It is expected that the results based on 
the velocity autocorrelation function are more accurate, because the evaluation at each fixed 
time s*, is based on an average over all times included in the simulation. In contrast, the 
mean square displacement is obtained from a single evaluation for each s*. This qualitative 
difference in accuracy has been confirmed by running multiple simulations for a few systems, 
and observing variations in the coefficient of self-diffusion as large as 10% when determined 
from the msd, whereas that determined from the vacf showed variations of less than 1%. 
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III. COMPARISON BETWEEN MD AND THEORY 



Figure |^ shows the ratio of the observed diffusion constants to the value predicted by the 
Enskog-level theory, as a function of a for densities n* = 0.1,0.25,0.5, and 0.75. For the 
lowest density, theory and simulation are in excellent agreement even at strong dissipation, 
as expected from earlier comparisons of simulation results and predictions of the Boltzmann- 
Enskog equation 0. As the density increases, deviations from the Enskog prediction are 
expected, since they are known to occur even at equilibrium {a = 1). In this respect, our 
results for a = 1 are quite consistent with previous studies of elastic hard spheres|0]. It is 
seen in the figure that deviations increase with decreasing a. For example, at n* = 0.25 
the deviations are less than 10% for a > 0.9, but increase rapidly to the order of 25% for 
a = 0.5. This behavior can also be observed in the normalized vacf itself which, in the 
Enskog approximation (here and below we refer to the first Sonine approximation to the 
Enskog theory), is given by exp(— cj^s*), with [|l| 

ul = ^-^^x{-T-f\ (14) 

where x is the pair correlation function for two particles at contact. In Fig. the logarithm of 
the normalized velocity autocorrelation function, C**{s*), is plotted as a function of uj^s*, for 
the case of elastic collisions, a = 1, and the several densities we have been considering. Here 
and below the choice has been made, somewhat arbitrarily but consistently, of truncating 
the vacf at lnC**(s*) ^ —4 in order to eliminate the tails which are dominated by noise. The 
data confirm findings of earlier studies that the Enskog theory gives a good description at 
low densities, but increasingly underestimates the diffusion constant, due to the neglect of 
correlated collisions and cage effects, at higher densities with a maximum deviation around 
n* = 0.5. Above this density, the neglected processes begin to cancel one another and near 
n* = 0.75, the diffusion constant crosses the Enskog prediction, and for higher densities is 
over-estimated by Enskog theory. Since the Enskog theory is exact at short times0], these 
effects appear as deviations from the simple exponential form of the vacf at longer times. 

Now consider a < 1. A similar plot. Fig. |^, for an inelastic system with a = 0.7, again 
shows good agreement with the Enskog prediction at the lowest density, but it exhibits 
much larger deviations than in the elastic case at higher densities, except at short times 
where the Enskog theory is exact, as already pointed out. Interestingly, the data clearly 
indicate a crossover to a slower constant decay rate at longer times. Of course, this behavior 
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increases the time integral of the velocity autocorrelation function, and it is responsible 
for the enhancement of the diffusion coefficient seen in Fig. ^ in contrast to the opposite 
behavior for high densities at a = 1. 

To characterize the deviation of the vacf from the simple exponential form, we start with 
an exact expression for the vacf based on the Zwanzig-Mori formalism (see Appendix |B]), 

d,,c:: is*) + uicz is*) + f ds' M{s* - s')c::{s') = o. (is) 

Jo 

where M{s) is known as the memory function. If M{s) is neglected the exponential decay of 
the above Enskog approximation is recovered, so the memory function incorporates all of the 
effects neglected in that approximation. A simple ansatz for this function as an exponential 
is qualitatively successful in modelling the vacf of fluids with elastic collisions (see p and 
references therein), and can be expected to work also for inelastic systems, when formulated 
in the dimensionless time s*. If we substitute M{s*) = M(0) exp(— As*) into Eq.(|T5D, and 
solve for the vacf with the boundary condition C**(0) = 1, the resulting model is 

C**{s*) = ~ ^ exp (-7+^iS*) + ^+ ~ ^ exp (-7-^1^*) , (16) 
7- - 7+ 7+ - 7- 

where the constants 7+ and 7_ can be related to M(0) and A. Figure |^ includes the result 

of fitting these two parameters to the data, being evident that the above model is able to 

capture the crossover from the Enskog behavior at short times to the slower relaxation for 

longer times. 

While the memory function model provides a framework for describing the results of the 
simulations, it does not explain them, since any effect not captured in the Enskog theory will 
give rise to a non-zero contribution to the memory function. In a previous study , similar 
deviations from the Enskog theory were found for the pressure of the system, and evidence 
was provided there suggesting that the large discrepancies between theory and simulation for 
dense, dissipative states may be due to the underlying hydrodynamic instability. Additional 
support for this possibility comes from Fig. ^ which shows the diffusion constants obtained 
from the simulation data for systems composed by 500 particles compared to those for 
systems composed by only 108 particles. For a density of n* = 0.25, there is a significant 
further increase of the diffusion constant in the large system relative to the small system 
as a decreases, with the diffusion constant of the 500 atom system growing to almost twice 
that of the 108 atom system at a = 0.5. At n* = 0.5, the enhancement is even larger, but 
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we know p[ that the system is in the unstable regime, at least for a < 0.7 and possibly 
for smaller values, so that most of the enhancement is undoubtedly due to the spontaneous 
formation of shear flow in the system. This suggests that in the smaller systems, although 
they are stable, there are present large fluctuations characteristic of the instability, such as 
spontaneous vortices that form and breakup. In this case, the impurity would find itself in 
a fluctuating local flow field which could enhance the velocity correlations. 

To test whether a local flow field plays a relevant role in the self-diffusion process, two 
different methods were used to calculate the autocorrelation of the impurity velocity relative 
to the instantaneous local flow field. In the first method, the instantaneous local flow field 
was calculated by dividing the simulation cell into 3^ = 27 cubic subcells. Then each atom's 
velocity relative to the instantaneous average velocity of the fluid in the subcell containing it, 
was used when calculating both the temperature and the vacf. For the 108 particle system, 
this procedure leads to a lowering of the measured value of the equilibrium diffusion constant 
by about 10%, which we believe to be a finite-size effect due to the removal of a significant 
fraction of the degrees of freedom of the system. For n* = 0.5 it also reduces the steady 
state temperature from nearly twice the initial temperature (see Fig. |I]) to about 1.5T**(0) 
thus demonstrating the size of the contribution of these fluctuations to the temperature. 
When the same procedure is applied to a system of 500 atoms, the shift in the equilibrium 
diffusion constant is negligible. The results obtained by this method are presented in Figs. 

In the graphs, the a-dependant diffusion constant has been scaled by its measured 
equilibrium value. The conclusion emerging from the flgures is that the a-dependence of 
the new diffusion constant is substantially closer to the Enskog form. However, this method 
may be criticized on the grounds that the statistics of the local velocity fleld are poor, since 
the calculations only involve a few atoms in each cell. We, therefore, also consider a second 
method which is specifically designed to eliminate only the part of the local flow due to the 
longest wavelength fluctuations in the system. A local flow fleld u is deflned by summing 
over the smallest Fourier components compatible with the size and shape of the simulation 
cell, 

u^{r)= J2 «4k)e^'^■^ (17) 

|k|>27r/L 

where is the volume of the cubic simulation cell. The Fourier components are determined 
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from the instantaneous velocities of the bulk fluid via 

«.(k) = ^$^t;.„e-'^■r^ (18) 

i 

This method is expected to give better statistics than the cell method discussed above. 
We then use this flow field to define the relative velocity used in the calculation of the 
temperature and the vacf. Figures |9|-|lT| also show the resulting diffusion constant when 
this method of subtracting the instantaneous fluctuations is used. From the figures follows 
that both methods are mutually consistent over the range 0.5 < a < 1, and that the 
agreement between the simulation results and the theoretical prediction from Enskog theory 
is significantly improved. The effects of subtracting the local flow field are much larger for 
the inelastic case due to the mechanism responsible for shear instability. 



IV. DISCUSSION 



In this paper, we have continued the discussion of diffusion in a model granular system 
begun in ref. The main motivation has been to show that transformation to the steady- 
state dynamics allows to carry over many standard methods of nonequilibrium statistical 
mechanics with relatively minor modifications. Here, the exact correspondence between the 
usual formulation of the dissipative hard-sphere model of granular fluids and the steady- 
state dynamics was exploited to formulate a particularly convenient simulation method 
which eliminates the need for additional complications such as exothermic boundary condi- 
tions. We have demonstrated that in the steady-state variables, the homogeneous cooling 
state exhibits standard diffusive behavior such as the linear increase of the mean-squared 
displacement with time, and the correspondence between the Einstein and Green-Kubo 
methods of determining the diffusion constant. In addition, the simplest approximation of 
Enskog kinetic theory provides excellent agreement at low densities for the whole range of 
inelasticity 0.5 < a < 1. However, although the formalism and concepts of the statistical 
mechanics of diffusion have been shown to give an adequate description of diffusion in the 
HCS, large quantitative deviations from the Enskog kinetic theory predictions have been 
observed even at relatively moderate densities. At the level of the vacf, it was shown that 
the deviations are principally due to a crossover from the Enskog time dependence, which 
is exact for short times, to another, slower decay which nevertheless also appears to be 
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exponential in form. The analytic form of the vacf was shown to be well approximated by 
modeling the memory function of the diffusive process by a simple exponential as has been 
used in early studies of memory effects in equilibrium fluids 0. 

Physically, we have shown by means of constrained simulations that the deviations from 
the Enskog model are primarily due to the effect of the longest wavelength velocity modes 
in the system. For sufficiently large systems, the shear mode is linearly unstable (in the 
steady-state picture) in the classic sense that the time constant associated with its decay 
goes to zero as the critical value of a is approached from above. A linear stability analysis 
in the steady-state variables!^ shows that the shear fluctuations of wavevector k decay 
exponentially with a time constant of 77* /c*^ — where if is the shear viscosity and the 
stars indicate quantities expressed in the reduced units of The dominant a-dependence 
of this expression comes from ~ (1 — a^) and the size of the system enters through the fact 
that the smallest non-zero value of the wavevector the system can sample is fcj^j^ = 27r/L* 
where L* is the longest dimension of the simulation cell. For fixed values of a < 1, there will 
always be a critical value of L* above which the system is unstable to shear fluctuations. 
For values of a above the critical value, the decay of these fluctuations is nevertheless slowed 
relative to equilibrium and we therefore reason that, even before the onset of the instability 
and even in systems too small to exhibit the instability, the slowing down of this mode means 
that long-lived, long-wavelength fluctuations are present. These fluctuations give rise to the 
observed enhancement of the diffusion constant. This effect is somewhat analogous to what 
is observed in turbulent systems since the removal of these modes in the calculation of the 
diffusion constant removed a large part of the discrepancy from the Enskog prediction. The 
size dependence of the deviations from the Enskog results supports this conclusion, since 
larger systems have a higher critical value of a because the wavevector can take on smaller 
values. We conclude that a complete theoretical description of diffusion in HCS will require 
a model for the memory function taking into account the slow relaxation of the shear modes 
by means, e.g., of a mode-coupling mechanism. 

V. ACKNOWLEDGMENTS 

The research of JWD was supported by National Science Foundation grant PHY 9722133. 
JJB acknowledges partial support from the Direccion General de Investigacion Cienti'fica 

14 



y Tecnica (Spain) through Grant No. PB98-1124. JL acknowledges support from the 
Universite Libre de Bruxelles. 



APPENDIX A: APPROACH TO STATIONARITY 

The average value of an observable A{r) for a general homogeneous state p{T,t) is given 
by Eq. (11) above which can be written 



{A; t) = J rfrp*(r )A({£q* (.) , t;(t)v* (.)}) (Al) 

and 

dsci: (s) = V* (s) , dsv* is) = ^CcsV* is) + L*v* (.) . (A2) 

The subscript on indicates explicitly that it is the dimensionless cooling rate associated 
with the HCS 

where i is the mean free path and Vhcsif) = 2Thcs{t) /tn is the HCS thermal velocity. The 
temperature Thcsit) obeys the equation 

dtThcsit) = -Chcs{Thcs{t))TrUt)- (A4) 
Now make the change of variables 

q-=q*, vr = (t/^7CL)v*, .* = (CL/^*)s, A**(.*)=A(.), (A5) 

to get 

{A-t) = I dr**p**{r*,s*)A{{iqr,v{t) (CL/^*) vD), (A6) 
with the definition dT** p** [T** , s*) = dr*p*{r*, s). In particular, the kinetic energy is 

{^mv'-t) = vLit) ia*Jw*)' I drp*{r*,s*)^mv**' 

= lrnvlMiCHj^*f{v**'-,s*y\ (A7) 
The corresponding temperature is 



2 1 /t* 



Tit) = ^{^mv';t) = 2TUt) ( ^ ) T**{s*) (A8) 



15 



where T**{s*) = d {v** , s*) . The time derivative of T**{s*) is then found to be 

(9,. - w*) T**{s*) = -w*^^T**{s*). (A9) 
This is still valid for a general homogeneous state. Now assume the existence of a scaling 



solution, which implies C{t) /Chcs{t) = \/T{t) /Thcs{t). Then 

{ds* - w*) T**{s*) = ~V2ChcsT**'^'\s*). (AlO) 
The solution is that given by Eq. (0) of the text. 

APPENDIX B: MEMORY FUNCTION MODEL 

The dimensionless velocity auto correlation function is defined by 

c:„(.)^(^(.)^r, = -!%=• (Bi) 

Using the notation in ref. this can be written in terms of the generators of the dynamics 
C:^{s) = [ dT* {e^'^^p) pL> = / rfF^e"^*^ (pL» • (B2) 



The detailed forms for the linear operators C* and C are given by Eqs. (36) and (42) of 
ref. JH but will not be required here. A projection operator is defined by 



PX = pl^^i, / dri^X. (B3) 



It follows then that 



PX{s)=pl,,i^CUs), (B4) 



with the choice 



X(.) = e-^^(pL.r)^). (B5) 
The equation of motion for X{s) is 

'a,+r)X(s) = (B6) 
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and a closed equation for PX{s) is obtained by operating on this equation with P and 
Q = 1 — P to get the pair of equations 

(ds + PC* P^ PX (s) = -PTQX{s), 

(d, + QTq) QX (s) = -QTPX {s) . (B7) 

Solving formally for QX (s) in the second equation and substituting into the first gives the 
desired closed equation for PX(s)[|l^ 

(^ds + PZ*P^ PX (s) - ds'e-^'^*^^'~'"^PTQTPX {s') = 0. (B8) 

Use of ( |B4| ) gives the corresponding equation for C*^(s) 

d^ClM + ^iClXs) + rds'M{s-s')C:,{s') = 0, (B9) 

^0 

with the definitions 

= {{C*^)^r (BIO) 

and 

M{s) = - j dri,Te-'^'^*'^'QTpl^^ij. (Bll) 

The definition for uji is the same as that for the Enskog approximation discussed in |]l| 
and approximated in eq. (|l^) above. Finally, it is clear that structurally identical expres- 
sions would have been obtained had we performed this derivation using the scaled variables 
introduced in Eq. (|^) of the text thus justifying the expression in Eq. (pIS]). 
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FIG. 1: Time evolution of the dimcnsionless steady temperature T** as a function of scaled di- 
mensionless time s* for n* = 0.5. The lower curve is for a = 1 , the middle for a = 0.7, and the 
upper for a = 0.5. 
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FIG. 2: Average mean square dimensionless displacement (msd) as a function of the dimensionless 
time in the steady-state dynamics for n* = 0.5. The lower curve is for a = 1 , the middle for 
a = 0.7, and the upper for a = 0.5. 
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FIG. 3: Normalized velocity autocorrelation function C{s*) as a function of the dimensionless 
scaled time s* in the steady-state dynamics for n* = 0.5. The values of the coefficient of restitution 
are a = 1 (lower curve), a = 0.7 (middle curve), and a = 0.5 (upper curve). 
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FIG. 4: Diffusion constant as a function of alpha for n* = 0.1 (circles) and n* = 0.5 (diamonds), as 
determined from the mean-squared-displacement (open symbols) and the velocity autocorrelation 
function (full symbols). The values are normalized with those obtained for a = 1 
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FIG. 5: Ratio of the diffusion constant determined by MD to the predicted value based on the 
Enskog theory, as a function of a, for n* = 0.1 (full circles), n* = 0.25 (full triangles), n* = 0.5 
(open circles), and n* = 0.75 (open triangles). The lines are a guide to the eye. 
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FIG. 6: Negative logarithm of the normahzed vacf C** as a function of time for equihbrium for a 
system in equihbrium (a = 1). Symbols are as in Fig. |5[ The lines are a guide to the eye, except 
the full line, which indicates the Enskog result. 
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FIG. 7: Negative logarithm of the normahzed vacf C** as a function of dimensionless time for the 
steady-state dynamics with a = 0.7. Symbols are as in Fig. The lines are the results of fits of 
the memory function model with an exponential kernel as discussed in the main text. The lowest 
density results are indistinguishable from the Enskog prediction on this scale. 
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FIG. 8: Diffusion constants for n* = 0.25 (circles) and n* = 0.50 (squares) as a function of a, for 
108 atoms (full symbols) and 500 atoms (open symbols). The lines are a guide to the eye. 
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FIG. 9: Diffusion constants for n* = 0.25 as a function of a for 108 atoms, calculated from the 
Green-Kubo relation (error bars and connecting line), by removing the longest Fourier modes 
(diamonds), and by using the cell-method of computing the local velocity field (circles). The heavy 
line is the prediction from Enskog theory. 
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FIG. 10: The same as Fig. | for n* = 0.50. 
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FIG. 11: The same as Fig. ^ for n* = 0.75. 
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